%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% (1.a) Forecasting future employment -> 2040
%%        using different age-specific ret-rates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   NB. we need cells for this: we want age-specific 
%        ages to take account of 
%               (i) diff rhos across ages
%               (ii) compositional changes over time of age structure
%
%       HERE DO NOT attempt
%               (iii) diff welfare costs across ages
%                   (this diff as suggests expectations about future rho..)
%
% Get transition and wage distribution parameters from Bestimcell
% Important to run Bestimcell *before* all else  => starts with "clear" command. 
% Note that Bestimcell contains age-specific rho-vector for cells 5,6,11,12
%


%% Which Sample & ROBUSTNESS TEST 

sample=2;
% Asimulage will only work with different age groups - this requires
% sample=2 also in Bestimcell - check there is you have problme
Bestimcell

% relative job search efficiency on-the-job with EGS/WIS 
% compared to unemployed post-coal job searchers relative-search-efficiency-J2J vs. U2J
% Default (robust1=0) is that with EGS/WIS, people find as many jobs as unemployed workers
% robust1=0    % benchmark:
% robust1=1;   % job search efficiency on adoption of WIS only 50% of standard
% robust1=2;   % if set robust1=2: firms reduce wages to minwage as result of WIS
 robust1=3;   % if set robust1=3: workers with twice job-finding efficiency in WIS
% robust1=4;   % if set robust1=4: workers with 1/4th job-finding efficiency in WIS  

if robust1==1
    disp("Robustness check: With wage insurance, j2j transitions: 50% of unemployed")
    setj2j=0.5;
elseif robust1==0
    setj2j=1;
    disp("Default without wage insurance, j2j transitions: 100% of unemployed")
elseif robust1==2
    setj2j=1; 
    disp("Robustness check for revise & resubmit:  benchmark search efficiency under WIS, minwage-job offers")
elseif robust1==3
    setj2j=2; 
    disp("Robustness check for revise & resubmit:  search efficiency twice as high as unemp")
elseif robust1==4
    setj2j=0.25; 
    disp("Robustness check for revise & resubmit:  search efficiency only one fourth of unemp")
end

%% Parameters
% Cell values of workforce in lignite in 2017 are reported in 
% results/${sample}/5_sample2_lig2017.xlsx
% (produced by section (7bii) of 5estimations.do)
% using 5_sample2_lig2017.xlsx, proportion of women is:
cell.propfem=1292/11625;

% years for which we want to make predictions
age.year=(2018:2040);
age.Yx=length(age.year);

% for WIS/EGS Policy simulation: job-to-job-with-egs-finding-rate is
% out-of-unemp-finding-rate-after-coal
cell.j2jd  =  setj2j.*param.lamNC/30.5;
cell.j2jm  =  setj2j.*param.lamNC;
cell.j2jy  =  setj2j.*(1-(1-cell.j2jm).^12);

% do not use values for women - rather use "female multiplier" to increase
% welfare costs calculated for men (small sample sizes create challenge)
cell.femult=1/(1-cell.propfem);

% number of cells
cell.Nx=length(param.wCmu);

% number of age groups
age.age=(18:65);
amin=18;
agenorm=17;

% for all groups apart from for ages 50-65, assume zero retirement
% for groups 5,6 and 11,12 (i.e. the cells with individuals aged 50-65)

% calculate yearly age-specific retirement rates
param.rhoy=zeros(65-agenorm,cell.Nx); % set size
param.rhoy(1:(49-agenorm),:)=0;
% these give rho-parameters for ages 50-65
param.rhoy(50-agenorm:65-agenorm,5) = [param.rhoy5];
param.rhoy(50-agenorm:65-agenorm,6) = [param.rhoy6];
param.rhoy(50-agenorm:65-agenorm,11)=[param.rhoy11];
param.rhoy(50-agenorm:65-agenorm,12)=[param.rhoy12];

% calculate daily age-specific retirement rates
param.rhod=zeros(65-agenorm,cell.Nx); % set size
param.rhod(1:(49-agenorm),:)=0;
% these give rho-parameters for ages 50-65
param.rhod(50-agenorm:65-agenorm,5) = [param.rhod5];
param.rhod(50-agenorm:65-agenorm,6) = [param.rhod6];
param.rhod(50-agenorm:65-agenorm,11)=[param.rhod11];
param.rhod(50-agenorm:65-agenorm,12)=[param.rhod12];

% (1) Rate of natural ageing. 
% initial(2017) age-distribution in the 4 different cells
% get this from 5agesN.m produced by 5estimations, called by Bestimcell.m 
age.Nx=[param.ageN1 param.ageN2 param.ageN3 param.ageN4 param.ageN5 param.ageN6 param.ageN7 param.ageN8 param.ageN9 param.ageN10 param.ageN11 param.ageN12];
% matrix of dimensions (ages30-65,cells5-6-11-12,years)
pop=zeros(size(age.Nx,1),size(age.Nx,2),age.Yx);
popj2j=zeros(size(age.Nx,1),size(age.Nx,2),age.Yx); % preallocate
pop(:,:,1)=age.Nx;
basepop=zeros(size(pop));
agedim=18:65;

%% Now calculate simplified version of Entgeltsicherung-effect:
% assume all jobs are accepted... increase outflow by lambda_NC
% (see Apolicy for calculation of associated costs)

setegs = 1;
if setegs==0
    disp("No-WIS setting: projecting future workforce sizes without WIS")
elseif setegs==1
    disp("WIS setting: projecting future workforce sizes with & without WIS")
end

% always run first without egs - later add run with egs=1 if setegs==1
egs.run = 0;
while egs.run < setegs + 1
if egs.run==1
    disp("Press enter to now project future workforce sizes including EGS")
    pause
end

% We now first define rates of natural ageing that will shift
% individuals between different age groups. 

% We are particularly interested in the four cells (5,6,11,12) that have over-50s
% but we also need to take into account inflow between cells:
% People from cells (1,2) age to become (3,4) and (3,4) become (5,6)
% Similarly, from cells (7,8) ppl move to (9,10) and from (9,10) to (11,12)
 % preallocate age- & cell-specific rates of leaving job
outrated=NaN(length(agedim),12);
outratey=NaN(length(agedim),12);  
outrate=NaN(length(agedim),12);  
j2jtrans=NaN(length(agedim),12);  
j2jd=NaN(length(agedim),12);  
j2jy=NaN(length(agedim),12);  
for c=1:12
    % for different projection years in future
    for y=2:2040-2017
        year=2017+y;
        %    population is matrix (age,cell,year)
        %    startpop is before outflow & inflow
        % no young 18-teen years olds joining: workforce ageing...
        basepop(i,c,y) = pop(i,c,y-1);
        for i=2:length(agedim)
            %   ageflow (1/3) simply shifts all pop one age-year with calendar-year on 
            basepop(i,c,y) = pop(i-1,c,y-1);
            %   ageflow (2/3) shifts 30-year olds across cells
            basepop(31-agenorm,3,y)=pop(30-agenorm,1,y-1);% shift 30-year olds from 1=>3 
            basepop(31-agenorm,4,y)=pop(30-agenorm,2,y-1);% shift 30-year olds from 2=>4
            basepop(31-agenorm,9,y)=pop(30-agenorm,7,y-1);% shift 30-year olds from 7=>9
            basepop(31-agenorm,10,y)=pop(30-agenorm,8,y-1);% shift 30-year olds from 8=>10
            %   ageflow (3/3) shifts 50-year olds across cells
            basepop(50-agenorm,5,y)=pop(49-agenorm,3,y-1);% shift 50-year olds from 1=>3 
            basepop(50-agenorm,6,y)=pop(49-agenorm,4,y-1);% shift 50-year olds from 2=>4
            basepop(50-agenorm,11,y)=pop(49-agenorm,9,y-1);% shift 50-year olds from 7=>9
            basepop(50-agenorm,12,y)=pop(49-agenorm,10,y-1);% shift 50-year olds from 8=>10
            % clean age categories of false natural ageing:
            % in 18-30 y.o. cells, no older ppl
            basepop(31-agenorm:end,1:2,y)=0; 
            basepop(31-agenorm:end,7:8,y)=0; 
            % in 31-49 y.o. cells, no younger/older ppl
            basepop(1:30-agenorm,3:4,y)=0; 
            basepop(1:30-agenorm,9:10,y)=0; 
            basepop(50-agenorm:end,3:4,y)=0; 
            basepop(50-agenorm:end,9:10,y)=0; 
            % in 50+ y.o. cells, no younger ppl
            % NB. if we use sum to calculate workers, this also implies
            % outrate=100% for 65-year olds...
            basepop(1:49-agenorm,5:6,y)=0;
            basepop(1:49-agenorm,11:12,y)=0;
            if egs.run==0
                outrate(i,c)=param.rhoy(i,c);
            elseif egs.run==1
            % here using daily rate to evaluate the ES/WIS policy to increase precision 
            % similar results with: outrate(i,c)=min(param.rhoy(i,c)+cell.j2jy(c),1); 
            % (note we have in cell.j2j also included the (modulable) search-efficiency
            % term setj2j (comparing search effort with EGS/WIS to
            % unemployed job-search).
            % here only adding job-to-job mobility (and entitlement to the
            % programme) for individuals up to age of 60
                if i < 60-agenorm
                    outrated(i,c)=param.rhod(i,c)+cell.j2jd(c);
                    j2jd(i,c)=cell.j2jd(c);
                elseif i>=60 - agenorm
                    j2jd(i,c)=0;
                    outrated(i,c)=param.rhod(i,c);
                end
                % generating yearly transition rates from the daily values
                j2jy(i,c)=1-(1-j2jd(i,c)).^365;
                outratey(i,c)=1-(1-outrated(i,c)).^365;
                outrate(i,c)=min(outratey(i,c),1); 
            end
            % now use egs.run-specific outrate to calculate 
            % change in population across different ages & years
            pop(i,c,y) = basepop(i,c,y).*(1-outrate(i,c));
            % specifically for the egs=1 runs, also count job-to-job movers
            if egs.run==1
            popj2j(i,c,y) = basepop(i,c,y).*(j2jy(i,c));            
            end
        end
    end
end

% population is matrix (age,cell,year)
% create matrix of cells x years (summing all ages)
% do this for (i) summing across cells for total workforce projection over simulation time
%             (ii) identify number of J2J movers every year as difference
%                  in workforces across EGS and non-EGS simulations.
if egs.run==0
    cell.pa=round(squeeze(sum(pop,1)));
    elseif egs.run==1
    cell.paegs=round(squeeze(sum(pop,1)));
    popj2j=round(squeeze(sum(popj2j,1)));
% see below: will create vector of sum of all cells- total workforce across years to project workforce size:
end

% welfare costs of coal exit in diff groups
cell.wfcpp=res.wfcpp;

% How to treat macro conditions? Assume good macro
% conditions into the future - pop in bad macro 
% is zero for year 2017 across all regions already anyway.
%
% welfare cost in every year 
% - applying pop-average welfare costs to pop declining with pop-mean rho
if egs.run==0
    cell.wfc=cell.femult.*cell.wfcpp'*ones(1,age.Yx).*cell.pa;
    cell.wfctot=sum(cell.wfc,1)/1000000;
elseif egs.run==1
    cell.wfc=cell.femult.*cell.wfcpp'*ones(1,age.Yx).*cell.paegs;
    cell.wfctotegs=sum(cell.wfc,1)/1000000;
end

disp("welfare costs are increased by female multiplier, which increases costs by (pct)")
disp((cell.femult-1)*100)

% Underlying Assumption: Welfare costs for women are average of those for men.
if egs.run==0
    disp("Evolution in workforce over time (no WIS/EGS)")
    cell.tot=sum(cell.pa,1); 
    disp("Projected workforce by 2025")
    disp(cell.tot(8)); pop25base=cell.tot(8);
    disp("Projected workforce by 2030")
    disp(cell.tot(13)); pop30base=cell.tot(13);
    disp("Projected workforce by 2038")
    disp(cell.tot(21)); pop38base=cell.tot(21);
    disp("Projected workforce by 2040")
    disp(cell.tot(23)); pop40base=cell.tot(23);
elseif egs.run==1
    disp("Evolution in workforce over time with WIS/EGS policy")
    cell.totegs=sum(cell.paegs,1);
    disp("Projected workforce by 2025")
    disp(cell.totegs(8)); pop25egs=cell.totegs(8);
    disp("Projected workforce by 2030")
    disp(cell.totegs(13)); pop30egs=cell.totegs(13);
    disp("Projected workforce by 2038")
    disp(cell.totegs(21)); pop38egs=cell.totegs(21);
    disp("Projected workforce by 2040")
    disp(cell.totegs(23)); pop40egs=cell.totegs(23);
end

disp("Welfare costs (m�) of exit in year 2017 based on cells & female multiplier")
disp(sum(cell.wfc(:,1)/1000000)); 
disp("Welfare costs (m�) of exit in year 2030 based on cells & female multiplier")
disp(sum(cell.wfc(:,13)/1000000)); 
disp("Welfare costs (m�) of exit in year 2038 based on cells & female multiplier")
disp(sum(cell.wfc(:,21)/1000000)); 
disp("Welfare costs (m�) of exit in year 2040 based on cells & female multiplier")
disp(sum(cell.wfc(:,23)/1000000)); 

egs.run=egs.run + 1;
end
% ends with egs.run=2


%% (1) Influence of WIS (EGS) on workforce numbers

disp("Reduction in workforce due to EGS by 2025 (%)")
disp((pop25base-pop25egs)./pop25base)
disp("Projected reduction (%) in workforce due to EGS by 2030")
disp((pop30base-pop30egs)./pop30base)
disp("Projected reduction (%) in workforce due to EGS by 2038")
disp((pop38base-pop38egs)./pop38base)
disp("Projected reduction (%) in workforce due to EGS by 2040")
disp((pop40base-pop40egs)./pop40base)

%% Graphing the influence of WIS (EGS in filename) on workforce numbers
clf
a1=plot(age.year,cell.tot);L1='without WIS';
axis([2017 2040 0 12000])
xlabel('years')
ylabel('projected total coal workforce')
if robust1==0
    saveas(gcf,'../../../paper/figures/projworkforce','png');
elseif robust1==1
     saveas(gcf,'../../../paper/figures/projworkforrob','png');
end
    
if setegs==1
    hold on
    a2=plot(age.year,cell.totegs); L2='with WIS';
    legend([a1,a2],L1,L2)
    if robust1==0
    saveas(gcf,'../../../paper/figures/projworkforceegs','png');
    elseif robust1==1
    saveas(gcf,'../../../paper/figures/projworkforceegsrob','png');
    end    
end

%% (2.a) Calculating the influence of EGS on welfare costs

disp("Reduction in welfare costs due to EGS by 2025 (%)")
disp((cell.wfctot(:,8)-cell.wfctotegs(:,8))./cell.wfctot(:,8))
disp("Projected reduction (%) in welfare costs due to EGS by 2030")
disp((cell.wfctot(:,13)-cell.wfctotegs(:,13))./cell.wfctot(:,13))
disp("Projected reduction (%) in welfare costs due to EGS by 2038")
disp((cell.wfctot(:,21)-cell.wfctotegs(:,21))./cell.wfctot(:,21))
disp("Projected reduction (%) in welfare costs due to EGS by 2040")
disp((cell.wfctot(:,23)-cell.wfctotegs(:,23))./cell.wfctot(:,23))

disp("Absolute numbers: Projected reduction in welfare costs due to EGS by 2030")
disp((cell.wfctot(:,13)-cell.wfctotegs(:,13)))
disp("Absolute numbers: Projected reduction in welfare costs due to EGS by 2040")
disp((cell.wfctot(:,23)-cell.wfctotegs(:,23)))

%% (2.b) Graphing the influence of WIS=EGS on welfare costs

clf
a1=plot(age.year,cell.wfctot);
L1='welfare cost without WIS';
axis([2017 2040 0 5000])
xlabel('years')
ylabel('projected total welfare costs (m�)')
if (robust1==0 || robust1==2)
    saveas(gcf,'../../../paper/figures/projwfc','png');
elseif robust1==1
        % Robustness check with limited job search efficiency on-the-job.
    saveas(gcf,'../../../paper/figures/projwfcrob','png');
end
if setegs==1 
    hold on
    a2=plot(age.year,cell.wfctotegs); L2 = 'welfare cost with WIS';
    legend([a1,a2],L1,L2)
    if (robust1==0 || robust1==2)
    saveas(gcf,'../../../paper/figures/projwfcegs','png');
    elseif robust1==1
    % Robustness check with limited job search efficiency on-the-job.
    saveas(gcf,'../../../paper/figures/projwfcegsrob','png');
    end
end

%% How much does the EGS cost?

% (1) How high are wage differences - when they do occur.
% Group-specific LLH of difference of two log normally distributed vars
% Here: simulate for simplicity and draw from this distribution.

rng(77);
sim.Nx1=1000000;
sim.diff=NaN(1,sim.Nx1);wwmu=NaN(cell.Nx,1); wwsig=NaN(cell.Nx,1); wwneg=NaN(cell.Nx,1); 
for c=1:cell.Nx % simulate for all groups 
sim.wNC(c,:) = lognrnd(wNCmu(c),wNCsig(c),[sim.Nx1,1]);
if robust1==2
% test alternative assumption of non-coal wages being reduced by firms to 
% minimum wage
    minwage=398; % monthly minimum wage for full-time job in 2017
    sim.wNC=ones(cell.Nx,sim.Nx1).*minwage;
end
sim.wC(c,:) =  lognrnd(wCmu(c),wCsig(c),[sim.Nx1,1]);
sim.diff = sim.wC(c,:)-sim.wNC(c,:);
sim.diffplus = sim.diff(sim.diff>0);
% fit lognormal distribution to approximate difference of wNC-wC
wwpara=lognfit(sim.diffplus);
wwmu(c)=wwpara(1);
wwsig(c)=wwpara(2);
wwneg(c)=sum(sim.diff<0)/length(sim.diff);
% calculate proportion of time that [ wNC > wC ]; % Alternat: simulate & draw
% nopay=logncdf(0,wwmu,wwsig); % how much probability mass in negative region?
end

% (2) When do these costs occur: 
% How many people leave due to EGS policy? 
% Rate of job-finding: population of groups times group-specific J2J-finding rate
% In our simulation, j2j mobility is only part of outrate in EGS policy
% simulation, thus number of pepole leaving is:

% Create matrix (cell-N,simulationyear-N) with number of J2J-movers in every year
% job-moving probability as used as one of the components of outrate for egs=1

% (3) Combine (1) and (2) to get approximate costs:
% 5 years, here no departures out of the wage-insurance scheme
% (Could reduce costs by assuming that some workers move to other jobs
%  => this would create discount rate - e.g.(1-delta-lamda-rho))

% start todo delme 22 Jan 2021
months=5*12;
% (1) establish cell-specific costs per month and worker
for c=1:cell.Nx
    % for each cell, we have mean cost times number of individuals in that
    % year. (Mean consists of fraction "sim.cheapct" of individuals with 
    % zero costs (for who wNC>wC)
    egs.meanww(c)=exp(wwmu(c)+wwsig(c)/2); % mean wage difference(C vs NC)

    % (2) apply cell-specific costs p.m. & p.w. 
%     to changing simulated pop with positive wloss in different cells and years 
    for y=2:length(age.year)
        
    % egs.j2jN is matrix (c,y) of simulated population of movers 
    % now take into account that fraction have no costs (as wNC-offer > wC)
    egs.j2jplus(c,y)=round(wwneg(c).*popj2j(c,y)); % how many people have no wloss 
    egs.j2jloss(c,y)=round((1-wwneg(c)).*popj2j(c,y));
    % in every year and cell, total costs given by: 
    % costs per-person-per-month (wage difference) times persons times months
    egs.costot(c,y)=egs.meanww(c).*egs.j2jloss(c,y).*months;
    end
end

disp("Mean costs of WIS per entitled person")
%disp(months.*egs.meanww)
disp(mean(months.*egs.meanww))

disp("Total of ppl moving job to job")
disp(sum(sum(popj2j)))

disp("Fraction of J2J movers entitled to WIS per year")
disp(1-mean(wwneg))

disp("Number of ppl entitled to WIS")
disp(sum(sum(egs.j2jloss,1)))

egs.costpa=sum(egs.costot,1).*cell.femult;
disp("Costs of WIS per year in m� - using female multiplier")
disp(egs.costpa/1000000)

disp("total costs of WIS policy - using female multiplier")
disp(sum(egs.costpa/1000000))

disp("total costs of WIS (w.fem multi) coal exit 2030")
disp(sum(egs.costpa(1:13)/1000000))

disp("total costs of WIS (w.fem multi) coal exit 2038")
disp(sum(egs.costpa(1:22)/1000000))

disp("total costs of WIS (w.fem multi) coal exit 2040")
disp(sum(egs.costpa/1000000))

%% Output table on key EGS results
clabmategs = {'WIS costs per entitled person \& month', 'months entitled', 'projected WIS users','costs (m Euro) of EGS - coal exit 2030','costs (m Euro) of EGS - coal exit 2040'};
rlabmategs = {'Whole sample'}; % results by cell group in an appendix table
mategsbase(1) = mean(egs.meanww);
mategsbase(2) = months;
mategsbase(3) = sum((sum(egs.j2jloss,2)));
mategsbase(4) = sum(egs.costpa(1:13)/10^6,2);
mategsbase(5) = sum(egs.costpa(1:23)/10^6,2);
% output results to latex
path= '../../../paper/' ;
if robust1==0
    file='egscostbase.tex'; filename= [path file];
elseif robust1==1
    file='egscostbaserob.tex'; filename= [path file];
elseif robust1==2
    file='egscostbaserob2.tex'; filename= [path file];    
end
matrix2latex(mategsbase, filename, 'rowLabels', rlabmategs, 'columnLabels', clabmategs, 'alignment', 'c', 'format', '%-6.2f', 'size', 'small');

%% Output table EGS results for all cells
clabmategs = {'Monthly WIS costs p.p.', 'projected WIS users','WIS Costs: exit 2030 (m\euro\ )','WIS Costs: exit 2040 (m\euro\ )'};
rlabmatcells = {'\textbf{1 Low Edu, 18-30, high unem}','\textbf{2 Low Edu, 18-30, low unem}','\textbf{3 Low Edu, 31-49, high unem} ','\textbf{4 Low Edu, 31-49, low unem}','\textbf{5 Low Edu, 50+, high unem}','\textbf{6 Low Edu, 50+, low unem}','\textbf{7 Hi Edu, 18-30, high unem}','\textbf{8 Hi Edu, 18-30, low unem}','\textbf{9 Hi Edu, 31-49, high unem}','\textbf{10 Hi Edu, 31-49, low unem}','\textbf{11 Hi Edu, 50+, high unem}','\textbf{12 Hi Edu, 50+, low unem}'};
% for appendix table on results by group
mategscells=NaN(Nx,4);
mategscells(:,1) = egs.meanww;
mategscells(:,2) = sum(egs.j2jloss,2);
mategscells(:,3) = sum(egs.costot(:,1:13)/10^6,2); % coal exit 2030
mategscells(:,4) = sum(egs.costot(:,1:23)/10^6,2);  % coal exit 2040

if robust1==1
    disp("Robustness check WIS modelling: j2j mobility at 50% of unemployed")
    setj2j=0.5;
elseif robust1==0
    setj2j=1;
    disp("Default WIS modelling: j2j mobility at 100% of unemployed")
end

% output results to latex
path= '../../../paper/' ;
if robust1==0
    file='egscostcells.tex';
elseif robust1==1
    file='egscostcellsrob.tex'; 
elseif robust1==2
    file='egscostcellsrob2.tex'; 
end
filename= [path file];
matrix2latex(mategscells, filename, 'rowLabels', rlabmatcells, 'columnLabels', clabmategs, 'alignment', 'c', 'format', '%-6.2f', 'size', 'small');
% use this to get percentatge of young people getting help with EGS

    